home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_Kn.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.6 KB  |  241 lines
  1. /* specfunc/bessel_Kn.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_gamma.h>
  26. #include <gsl/gsl_sf_psi.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28.  
  29. #include "error.h"
  30.  
  31. #include "bessel.h"
  32.  
  33. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34.  
  35. /* [Abramowitz+Stegun, 9.6.11]
  36.  * assumes n >= 1
  37.  */
  38. static
  39. int
  40. bessel_Kn_scaled_small_x(const int n, const double x, gsl_sf_result * result)
  41. {
  42.   int k;
  43.   double y = 0.25 * x * x;
  44.   double ln_x_2 = log(0.5*x);
  45.   double ex = exp(x);
  46.   gsl_sf_result ln_nm1_fact;
  47.   double k_term;
  48.   double term1, sum1, ln_pre1;
  49.   double term2, sum2, pre2;
  50.  
  51.   gsl_sf_lnfact_e((unsigned int)(n-1), &ln_nm1_fact);
  52.  
  53.   ln_pre1 = -n*ln_x_2 + ln_nm1_fact.val;
  54.   if(ln_pre1 > GSL_LOG_DBL_MAX - 3.0) GSL_ERROR ("error", GSL_EOVRFLW);
  55.  
  56.   sum1 = 1.0;
  57.   k_term = 1.0;
  58.   for(k=1; k<=n-1; k++) {
  59.     k_term *= -y/(k * (n-k));
  60.     sum1 += k_term;
  61.   }
  62.   term1 = 0.5 * exp(ln_pre1) * sum1;
  63.  
  64.   pre2 = 0.5 * exp(n*ln_x_2);
  65.   if(pre2 > 0.0) {
  66.     const int KMAX = 20;
  67.     gsl_sf_result psi_n;
  68.     gsl_sf_result npk_fact;
  69.     double yk = 1.0;
  70.     double k_fact  = 1.0;
  71.     double psi_kp1 = -M_EULER;
  72.     double psi_npkp1;
  73.     gsl_sf_psi_int_e(n, &psi_n);
  74.     gsl_sf_fact_e((unsigned int)n, &npk_fact);
  75.     psi_npkp1 = psi_n.val + 1.0/n;
  76.     sum2 = (psi_kp1 + psi_npkp1 - 2.0*ln_x_2)/npk_fact.val;
  77.     for(k=1; k<KMAX; k++) {
  78.       psi_kp1   += 1.0/k;
  79.       psi_npkp1 += 1.0/(n+k);
  80.       k_fact    *= k;
  81.       npk_fact.val *= n+k;
  82.       yk *= y;
  83.       k_term = yk*(psi_kp1 + psi_npkp1 - 2.0*ln_x_2)/(k_fact*npk_fact.val);
  84.       sum2 += k_term;
  85.     }
  86.     term2 = ( GSL_IS_ODD(n) ? -1.0 : 1.0 ) * pre2 * sum2;
  87.   }
  88.   else {
  89.     term2 = 0.0;
  90.   }
  91.  
  92.   result->val  = ex * (term1 + term2);
  93.   result->err  = ex * GSL_DBL_EPSILON * (fabs(ln_pre1)*fabs(term1) + fabs(term2));
  94.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  95.  
  96.   return GSL_SUCCESS;
  97. }
  98.  
  99.  
  100. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  101.  
  102. int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result)
  103. {
  104.   n = abs(n); /* K(-n, z) = K(n, z) */
  105.  
  106.   /* CHECK_POINTER(result) */
  107.  
  108.   if(x <= 0.0) {
  109.     DOMAIN_ERROR(result);
  110.   }
  111.   else if(n == 0) {
  112.     return gsl_sf_bessel_K0_scaled_e(x, result);
  113.   }
  114.   else if(n == 1) {
  115.     return gsl_sf_bessel_K1_scaled_e(x, result);
  116.   }
  117.   else if(x <= 5.0) {
  118.     return bessel_Kn_scaled_small_x(n, x, result);
  119.   }
  120.   else if(GSL_ROOT3_DBL_EPSILON * x > 0.25 * (n*n + 1)) {
  121.     return gsl_sf_bessel_Knu_scaled_asympx_e((double)n, x, result);
  122.   }
  123.   else if(GSL_MIN(0.29/(n*n), 0.5/(n*n + x*x)) < GSL_ROOT3_DBL_EPSILON) {
  124.     return gsl_sf_bessel_Knu_scaled_asymp_unif_e((double)n, x, result);
  125.   }
  126.   else {
  127.     /* Upward recurrence. [Gradshteyn + Ryzhik, 8.471.1] */
  128.     double two_over_x = 2.0/x;
  129.     gsl_sf_result r_b_jm1;
  130.     gsl_sf_result r_b_j;
  131.     int stat_0 = gsl_sf_bessel_K0_scaled_e(x, &r_b_jm1);
  132.     int stat_1 = gsl_sf_bessel_K1_scaled_e(x, &r_b_j);
  133.     double b_jm1 = r_b_jm1.val;
  134.     double b_j   = r_b_j.val;
  135.     double b_jp1;
  136.     int j;
  137.  
  138.     for(j=1; j<n; j++) {
  139.       b_jp1 = b_jm1 + j * two_over_x * b_j;
  140.       b_jm1 = b_j;
  141.       b_j   = b_jp1; 
  142.     } 
  143.     
  144.     result->val  = b_j;
  145.     result->err  = n * (fabs(b_j) * (fabs(r_b_jm1.err/r_b_jm1.val) + fabs(r_b_j.err/r_b_j.val)));
  146.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  147.  
  148.     return GSL_ERROR_SELECT_2(stat_0, stat_1);
  149.   }
  150. }
  151.  
  152.  
  153. int gsl_sf_bessel_Kn_e(const int n, const double x, gsl_sf_result * result)
  154. {
  155.   const int status = gsl_sf_bessel_Kn_scaled_e(n, x, result);
  156.   const double ex = exp(-x);
  157.   result->val *= ex;
  158.   result->err *= ex;
  159.   result->err += x * GSL_DBL_EPSILON * fabs(result->val);
  160.   return status;
  161. }
  162.  
  163.  
  164. int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array)
  165. {
  166.   /* CHECK_POINTER(result_array) */
  167.  
  168.   if(nmin < 0 || nmax < nmin || x <= 0.0) {
  169.     int j;
  170.     for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;
  171.     GSL_ERROR ("domain error", GSL_EDOM);
  172.   }
  173.   else if(nmax == 0) {
  174.     gsl_sf_result b;
  175.     int stat = gsl_sf_bessel_K0_scaled_e(x, &b);
  176.     result_array[0] = b.val;
  177.     return stat;
  178.   }
  179.   else {
  180.     double two_over_x = 2.0/x;
  181.     gsl_sf_result r_Knm1;
  182.     gsl_sf_result r_Kn;
  183.     int stat_0 = gsl_sf_bessel_Kn_scaled_e(nmin,   x, &r_Knm1);
  184.     int stat_1 = gsl_sf_bessel_Kn_scaled_e(nmin+1, x, &r_Kn);
  185.     int stat = GSL_ERROR_SELECT_2(stat_0, stat_1);
  186.     double Knp1;
  187.     double Kn   = r_Kn.val;
  188.     double Knm1 = r_Knm1.val;
  189.     int n;
  190.  
  191.     for(n=nmin+1; n<=nmax+1; n++) {
  192.       if(Knm1 < GSL_DBL_MAX) {
  193.         result_array[n-1-nmin] = Knm1;
  194.         Knp1 = Knm1 + n * two_over_x * Kn;
  195.         Knm1 = Kn;
  196.         Kn   = Knp1;
  197.       }
  198.       else {
  199.         /* Overflow. Set the rest of the elements to
  200.      * zero and bug out.
  201.      * FIXME: Note: this relies on the convention
  202.      * that the test x < DBL_MIN fails for x not
  203.      * a number. This may be only an IEEE convention,
  204.      * so the portability is unclear.
  205.      */
  206.         int j;
  207.     for(j=n; j<=nmax+1; j++) result_array[j-1-nmin] = 0.0;
  208.         GSL_ERROR ("overflow", GSL_EOVRFLW);
  209.       }
  210.     }
  211.  
  212.     return stat;
  213.   }
  214. }
  215.  
  216.  
  217. int
  218. gsl_sf_bessel_Kn_array(const int nmin, const int nmax, const double x, double * result_array)
  219. {
  220.   int status = gsl_sf_bessel_Kn_scaled_array(nmin, nmax, x, result_array);
  221.   double ex = exp(-x);
  222.   int i;
  223.   for(i=0; i<=nmax-nmin; i++) result_array[i] *= ex;
  224.   return status;
  225. }
  226.  
  227.  
  228. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  229.  
  230. #include "eval.h"
  231.  
  232. double gsl_sf_bessel_Kn_scaled(const int n, const double x)
  233. {
  234.   EVAL_RESULT(gsl_sf_bessel_Kn_scaled_e(n, x, &result));
  235. }
  236.  
  237. double gsl_sf_bessel_Kn(const int n, const double x)
  238. {
  239.   EVAL_RESULT(gsl_sf_bessel_Kn_e(n, x, &result));
  240. }
  241.